Stepping stones to extirpation: Puma patch occupancy thresholds in an urban‐wildland matrix

Abstract Habitat loss and fragmentation are the leading causes of species range contraction and extirpation, worldwide. Factors that predict sensitivity to fragmentation include high trophic level, large body size, and extensive spatial requirements. Pumas (Puma concolor) exemplify these qualities, making them particularly susceptible to fragmentation and subsequent reductions in demographic connectivity. The chaparral‐dominated ecosystems surrounding the greater San Francisco Bay Area encompass over 10,000 km2 of suitable puma habitat, but inland waterways, croplands, urban land uses, and extensive transportation infrastructure have resulted in widespread habitat fragmentation. Pumas in this region now exist as a metapopulation marked by loss of genetic diversity, collisions with vehicles, and extensive human–puma conflict. Given these trends, we conducted a photo survey from 2017 to 2021 across 19 patches of predicted habitat and compiled a dataset of >6584 puma images. We used a logistic regression analytical framework to evaluate the hypothesis that puma patch occupancy would exhibit a threshold response explained by patch size, isolation, and habitat quality. Contrary to predictions, only variables related to patch size demonstrated any power to explain occupancy. On average, occupied patches were 18× larger than those where they were not detected (825 ± 1238 vs. 46 ± 101 km2). Although we observed pumas in patches as small as 1 km2, logistic regression models indicated a threshold occupancy probability between 300 and 400 km2, which is remarkably close to the mean male puma home range size in coastal California (~381 km2). Puma populations dependent on habitats below this value may be susceptible to inbreeding depression and human–wildlife conflict, and therefore vulnerable to extirpation. For species conservation, we suggest conflicts might be ameliorated by identifying the largest, isolated patches for public education campaigns with respect to management of domestic animals, and remaining connective parcels be identified, mapped, and prioritized for targeted mitigation.


| INTRODUC TI ON
Loss of demographic connectivity is a leading cause of species range contractions and extirpations, worldwide (Crooks et al., 2017;Jacobson et al., 2019;Semper-Pascual et al., 2021). Expansion of agriculture, urban land uses, and transportation infrastructure results in fragmentation, systematically reducing patch size, and increasing isolation and edge-area ratios. Anthropogenic barriers such as dams, fencing, or highways, can reduce or eliminate demographic connectivity in aquatic and terrestrial systems (Marschall et al., 2011;Seidler et al., 2014). High edge-area ratios in occupied habitat can simultaneously facilitate the spread of invasive species (Haddad et al., 2015) and promote human-wildlife conflict (Woodroffe & Ginsberg, 1998). These patterns interact to reduce the suitability of remaining habitats or isolate them altogether, resulting in non-linear trends in the rate of habitat loss and therefore community composition (Wilson et al., 2016). Yet, despite the ubiquity and acceleration of this problem, few extirpation thresholds of habitat area have been estimated for species of conservation concern.
Factors predicting vulnerability to extirpation have been studied extensively. Life history traits, niche specialization, and trophic level all contribute (Davidson et al., 2009), but isolation and habitat area have proven strong predictors of population persistence (Crooks et al., 2011;González-Suárez & Revilla, 2014). The effects of these factors are scale-sensitive, with large-bodied, obligate carnivores most likely to be impacted by current land-use trends because of extensive spatial requirements and low population densities (Stoner et al., 2018). Therefore, species at high trophic levels come into contact with hard boundaries at greater frequencies than those with smaller home ranges and resource needs (Woodroffe & Ginsberg, 1998). As such, retaining or restoring demographic connectivity among sub-populations has become one of the prevailing themes in wildlife management, with major conservation efforts focused on carnivores and/or migratory ungulates (e.g., USDI, 2018).
Pumas (Puma concolor) are one of the most broadly distributed mammals in the western hemisphere, inhabiting a wide variety of climatic zones and land-use types, including near-urban environments (Blecha et al., 2016, Riley et al., 2021, Stoner et al., 2021 Figure 1). As large-bodied felids, pumas are strict carnivores that depend on the presence and abundance of ungulate prey to support them (Pierce & Bleich, 2003). They occupy the highest trophic level, making them numerically rare in comparison with similar-sized omnivores and herbivores in the same systems. Although buffered from extinction by their large geographic range and ecological tolerance (Culver et al., 2000), pumas exemplify many of the qualities that make a species vulnerable to extirpation at local scales (Davidson et al., 2009;Purvis et al., 2000;Stoner et al., 2018).
Scientific investigations of pumas have occurred throughout the species' North American range. Despite exceptional dispersal abilities (Hawley et al., 2016;Stoner et al., 2013), research indicates that when combined with natural habitat patchiness, anthropogenic features can amplify fragmentation, thereby constraining puma movements (Stoner et al., 2013), or isolating subpopulations (Benson et al., 2019;Maehr et al., 2002). Some of the most compelling examples of this come from California and Florida (e.g., Beier, 1995;Maehr et al., 2002). The Florida panther (Puma concolor coryi) exists as a single, relict population with no connectivity to other extant populations, whereas California represents a bellwether for the effects of habitat fragmentation on species demographic connectivity (Benson et al., 2019;Dellinger et al., 2020). Indeed, several California puma populations number fewer than 100 individuals (e.g., Beier, 1993;Benson et al., 2020), and thus occur as sub-populations within a metapopulation context (Beier, 1993;Sweanor et al., 2000).
There is a growing concern among state wildlife agencies over the conservation value of small habitat patches (Fahrig et al., 2022), but as yet, there are no estimates of threshold values of habitat area to proactively identify populations at risk. Dellinger et al. (2020) calculated that 8000-15,000 km 2 of connected habitat were required to maintain an effective population size of 50 adult pumas and therefore genetic integrity. These authors identified five populations in California that did not meet this threshold, prompting policies

T A X O N O M Y C L A S S I F I C A T I O N
Applied ecology F I G U R E 1 Pumas occupy wildlands adjacent to major urban areas, but also traverse developed landscapes where they are vulnerable to various forms of human-wildlife conflict (photos courtesy of Steve Winter and Andy Forward). designed to provide greater protections for those subpopulations vulnerable to inbreeding and subsequent declines or extirpation.
Based on these trends, we set out to test the hypothesis that puma patch occupancy, and by proxy, local range contractions, would exhibit a threshold response to habitat fragmentation (Crooks, 2002).
We predicted that habitat patch size, isolation, and quality would best explain trends in puma presence (sensu MacArthur & Wilson, 1967), and that this would be influenced by patterns in land use and landcover. Furthermore, we predicted that area thresholds would correlate with puma spatial requirements typical for local environmental conditions.

| Study area
To evaluate this hypothesis, we estimated puma patch occupancy across the nine counties that comprise the greater San Francisco Bay Area in northern California, USA. The region measures 18,152 km 2 , and includes Sonoma, Napa, Marin and Solano Counties in the north, Contra Costa and Alameda Counties in the east, and San Francisco, San Mateo, and Santa Clara Counties in the south (Figure 2). Based on U.S. Census Bureau data published in 2020, the population of these counties was approximately 7.5 million people at the time of the study (https://www.census. gov/). Land uses range from minimally disturbed wildlands to dense urban areas, with suburban, exurban, agricultural, industrial areas, and open spaces comprising the gradient between wilderness and urban landscapes. The climate is Mediterranean, defined by warm, dry summers (10-33°C) and cool, rainy winters (2-18°C; 500-1200 mm precipitation; https://www.uscli mated ata.com/). Local climate varies spatially as a function of elevation and distance from the ocean. The marine fog belt maintains cooler and more consistent temperatures along the coast as compared to inland sites at the same elevation and latitude. Dominant plant communities reflect this climatic regime, and vary from chaparral shrublands near the coast, which grade into mixed oak woodlands and grasslands further inland. The area exhibits a high degree of natural fragmentation stemming from waterways that comprise the San Francisco Bay, Suisun Bay, and the Sacramento Delta.
Anthropogenic land uses are draped over these drainage patterns and serve to both amplify fragmentation and attract wildlife (Coon et al., 2020).

| Study patches
To identify habitat patches for field sampling, we used predictions of suitable habitat from Coon et al. (2020). This model used 4 years of camera trap data to predict occupancy probabilities for every 1 km 2 of the study area during both wet (winter) and dry (summer) seasons. The model was built with positive occupancy estimates for both forest cover (including evergreen, mixed, and deciduous forests) and distance to roads. For the current analysis, we sampled study patches from the dry season model by selecting the 1 km 2 pixels with the highest 50% real and extrapolated occupancy probabilities, and then grouped all 1 km 2 pixels adjacent to at least one other top-50% patch. Predicted occupancy patches for pumas were larger in the dry season, and therefore we used this delineation to draw habitat fragment boundaries. The model produced occupancy predictions for 51 individual patches in the study area, ranging in size from 1 to 4000 km 2 . However, 32 of these patches were eliminated from the sample frame due to inaccessibility (private or military lands), or resource limitations which precluded sampling of the smallest, most isolated patches or those dominated by wetlands with little structural vegetative cover ( Figure 2).

| Puma occupancy
Our response variable was puma presence (or detection) within a given habitat patch, which we used as an index of occupancy.
We did not use formal occupancy models per se (sensu MacKenzie et al., 2003), but used the term in the vernacular, meaning occurring in a place. To measure this, we sampled photos from a master database of more than 329,000 images derived from 483 remote camera placements conducted from 2017 to 2021. Trail camera images were cataloged by project biologists and volunteers, and candidate puma photos were validated by the author (CACC).
From this dataset, we compiled 6584 images of pumas from across the study area. Within selected habitat patches, initial camera deployments were placed strategically to maximize the probability of detecting wildlife on public lands and sites where we had permission to access private lands. Our criteria included the following: cameras were always ≥2 km apart, using common trigger speeds and motion sensors to maintain consistency; and no scents or lures were used to repel or attract wildlife. Camera deployment length varied by site, and as such, there was no maximum length of deployment. Instead, we defined a minimum threshold for justifying puma absence. Camera stations were established at densities of approximately 1 camera per 5 km 2 (0.2 cameras per km 2 ), and were active for 1-1974 days (mean = 214, SD = 269 days). Within this dataset, 180 camera placements had one or more puma detections and 94 had at least two detections (mean number detections at cameras with ≥1 detection = 11.2, SD = 26.7). Puma occupancy of a patch was determined through collection of at least one piece of verifiable evidence. These data came primarily from the long-term camera survey conducted by our non-profit organization (Felidae Conservation Fund, Bay Area Puma Project), and were supplemented with images provided by colleagues at state, county, and non-profit agencies, photos with geotags submitted by email or to our organizational citizen scientist online database (www.BAPP.org/sight ings-map), and publicly available observations of pumas marked by other organizations coincident with our camera survey (2016-2020; Table 2). Detections were not treated as independent.
Presence/absence surveys are prone to false negatives, in which putative absences do not preclude the possibility that animals were present but went undetected during the sampling interval. To derive a rule set for determining absence, we used data from cameras that were in the field for at least 115 days and produced one or more puma detections. Preliminary calculations indicated 90% of confirmed detections occurred within the first 124 days of sampling. Additionally, at camera sites in which a puma was detected on at least two occasions, mean latency between detections was 80 days. We averaged these two values to set a minimum density and sampling duration to evaluate presumed absences. Specifically, we multiplied 0.2 cameras/km 2 by 102 days to get a minimum required 20.4 camera-days/km 2 in each patch.
Although false negatives were still possible, this criterion provided a consistent, repeatable standard for justifying puma absence from a given habitat patch.

| Predictor variables
To assess how patch size affected patterns in puma occupancy, we tested three variables describing different aspects of area, including patch size (km 2 ), perimeter (km), and the perimeterarea ratio. All were calculated using ArcGIS V. 10.7 software (Environmental Systems Research Institute). We also tested secondary variables hypothesized to affect patch occupancy. These included "isolation," defined here as the distance between patches weighted by patch size (MacGarigal & Marks, 1995; Table 1); and "naturalness"-a continuous index of habitat quality that uses landcover, housing density, road presence, and traffic volume to scale anthropogenic disturbance (Theobold et al., 2012). Isolation and naturalness were calculated for each patch and an accompanying 1 km buffer. The naturalness buffer was used as a proxy for ecotones, that is, forest-grassland edges, which provide foraging habitat for mule deer (Odocoileus hemionus), and stalking cover for pumas (Holmes & Laundre, 2006). We then calculated the percent F I G U R E 2 Outline of nine San Francisco Bay Area counties in California that constitute the study area. Color-coded habitat patches are based on (1) whether the patch was enrolled in the study, and (2) whether or not a puma was detected during the sampling period (2017)(2018)(2019)(2020)(2021).
of landcover types that influence puma occupancy and prey vulnerability, including grassland/agriculture, and forest cover (Coon et al., 2020). Lastly, we evaluated several variables that have particular relevance in a study area characterized by extensive habitat alterations, such as non-native vegetation and impermeable surfaces. These included land ownership (% private), percent freshwater, percent developed, and relative road length. All variable definitions, units, and sources are detailed in Table 1.
TA B L E 1 List of variables, measurement units, source data, transformations, and t-test results for factors hypothesized to explain puma patch occupancy (San Francisco Bay Area, CA, 2017

| Analyses
We used Welch's two-sample t-tests to assess the significance of individual predictor variables on puma detection. T-tests were run with the t.test() function from the Stats package in program R (R Core Team, 2013). We then built a global multiple logistic regression model using variables with significant and marginally significant t-test results, while accounting for correlation between variables.
Logistic regression models were analyzed with the generalized linear model (glm()), also from the Stats package, using a binomial distribution and a logit link, with puma presence (1) or absence (0)

| RE SULTS
Of the 51 candidate patches identified from the puma habitat model, we were able to conduct field surveys, or obtain photographic evidence from 19 patches. Eleven patches produced evidence of being used by at least one puma during the sampling interval ( Figure 2, of ownership within patches, it was not always possible to systematically sample entire patches, hence the reason we report camera numbers and densities within study patches (Appendix 1). At a camera density of 0.2 cameras/km 2 for 102 days, we required a minimum of 20.4 (camera-days)/km 2 to confirm the absence of pumas in the study patch, which was met for all 19 focal patches.
Variables that provided the strongest evidence for differences in occupancy were mean patch size, patch perimeter length, perimeterarea ratio, percent grassland/agriculture, and percent anthropogenic development (p < .15; Table 1, Figure 3). We found no statistical differences between occupied and unoccupied patches with respect to isolation, naturalness, naturalness within a 1 km buffer, percent forest cover, percent private ownership, percent fresh water, and summed length of roads (p > .30; Table 1, Figure 3).
We used a combination of t-test results and cases of multicollinearity to determine which variables to include in the global model.
The global model included: patch size, patch perimeter, percent grassland/agriculture, percent development, percent fresh water, and road length within a patch, but we were unable to use all variables with significant or marginally significant t-test results (p < .10) due to multicollinearity among predictors. For example, patch size, perimeter, and perimeter-area ratio all had significant t-tests ( were 18× larger than those lacking detections (825 ± 1238 km 2 vs. 46 ± 101 km 2 ). Based on the simple logistic regression model with log (patch size) as the sole predictor, we identified a decrease in detection probability at patch sizes below log values of 2.5, or ~300-400 km 2 (n = 19, est. = 0.43; z = 1.81; p = .07; Figure 4).

| DISCUSS ION
Our objective was to evaluate patterns in puma habitat patch occupancy in a region defined by extensive levels of natural and anthropogenic fragmentation. We defined occupancy broadly as presence within a patch under the assumption that detection of pumas connotes some conservation value, even if only used as a stepping stone within an array of larger patches (Beier, 1995;Lynch, 2019). As expected, variables related to patch size (area, perimeter length, and perimeter-area ratio) displayed consistently strong relationships with puma occupancy. Pumas were detected in patches as small as 1 km 2 , but logistic regression results indicated a threshold value of 300-400 km 2 , suggesting that patches below this size were unlikely to harbor pumas. In a similar analysis, Crooks (2002) reported that probability of puma detection was lowest in patches <15 km 2 , and highest in those exceeding 100 km 2 . The author made clear that population viability was questionable at the bottom end of those estimates. Landcover was largely uninformative, with pumas less likely to be detected in patches where buildings and pavement accounted for more than 50% of the area. However, this is likely an artifact of our sampling scheme, in which we used selection criteria that minimized variation in landcover characteristics, a priori. Patches with no puma detections were relatively small, sparsely vegetated, or dominated by urban land uses. Contrary to expectation, isolation and naturalness had no discernable effects on puma detection. This presents an apparent contradiction in the literature that is mirrored in our results: pumas prefer large, natural spaces and tend to avoid humans even in highly populated areas (Benson et al., 2016), yet they persist where prey resources are available, regardless of patch naturalness.
We delineated habitat patches based on the model presented in Coon et al. (2020), in which forest cover and road infrastructure were the strongest positive and negative predictors of habitat quality, respectively. At the scale of the individual patch neither of these variables were important in predicting puma detection. We suspect the lack of a relationship is related more to sampling criteria than to any inherent behavioral tendencies expressed by pumas. Patches targeted for sampling met some minimum values with respect to roads and forest cover, and as such, did not capture the full range of variation that exists across the greater study area. Moreover, the occupancy model used to create the patches in this study (Coon et al., 2020) was restricted to the 9-county Bay Area, and as such, some of the edge patches based on county boundaries may be larger than estimated here.
The only landcover variable that positively impacted occupancy was the percentage of open habitat, such as grassland, pasture, or cultivated land. Puma habitat models consistently identify open, flat, or sparsely vegetated habitats as underused relative to availability, presumably because these cover types are incompatible with their stalk-and-pounce hunting style (Dickson & Beier, 2002;Logan & Irwin, 1985;Smereka et al., 2020). Our results largely support this generality, yet these same cover types and associated edges are preferred by black-tailed deer (Odocoileus hemionus columbianus), the primary puma prey species in this region (Allen, 2014;Hopkins, 1989), as well as various synanthropic species (Bateman & Fleming, 2012). Pumas are successful at hunting in forest-grassland ecotones (Holmes & Laundre, 2006), and therefore the weak but positive correlation between puma presence and open habitat may be an optimization of these constraints.

TA B L E 3
Values of variable importance on a 0 to 1 scale calculated from ranked models from all possible models generated by a global multiple logistic regression or some subset based on ΔAICc from the best model. Note: Variables included patch perimeter, patch development, and percent pasture and agriculture (denoted with an * above). together, this suggests that this patch may already be sufficiently isolated to reduce immigration and may therefore be vulnerable to extirpation.

F I G U R E 4
Pumas have been documented traveling through residential and urban environments (Riley et al., 2021;Suraci et al., 2020), but there are no examples of them occupying these areas indefinitely (Beier et al., 2010). Thus, the question still remains as to how animals are moving among the more isolated patches, given that indices of isolation had no effect on puma detection. Two recent analyses may provide some insights to this question. Suraci et al. (2020) studied pumas in the south-western portion of the study region and suggested micro-scale movement decisions based largely on attraction to vegetative cover and avoidance of urban landcover types. Other models of mammalian navigation suggest that in areas of high relief, transient animals may survey areas within line-of-sight prior to making extensive dispersal or migratory movements (Berger et al., 2022;Sweanor et al., 2000). Taken Beier (1996) estimated that individual patches of 1000-2200 km 2 would secure viability of a subpopulation at multi-decadal scales. More recently, calculations by Dellinger et al. (2020) suggested that 10,000 km 2 of contiguous habitat would be required to maintain an effective population size of 50 adult pumas to mitigate the effects of inbreeding. Puma social organization is characterized by a resident male overlapping two to five often-related adult females (Logan & Sweanor, 2010). In the Mediterranean climates of coastal California puma home range size varies widely, but averages 153 and 381 km 2 for females and males, respectively (Allen, 2014;Dickson & Beier, 2002;Hopkins, 1989 Table 2, the thresholds we present here might be used as indices of extirpation vulnerability and for prioritizing crossing structures between patches with the greatest connective value (e.g., Burdett et al., 2010;Crooks et al., 2011). Second, small patch size and high edge-area ratios can result in frequent conflict and high mortality (Benson et al., 2023;Woodroffe & Ginsberg, 1998).

| MANAG EMENT IMPLIC ATIONS
To reduce the potential for conflict associated with domestic animal depredation, isolated, but occupied patches should be targeted for public outreach and education activities (Vickers et al., 2015).
Lastly, private lands are highly vulnerable to development but are critical for preserving the connectivity that still exists. To the extent possible, improved land-use planning and permanent protection of suitable, connective habitat (Zeller et al., 2017) should be identified, mapped, and prioritized for targeted conservation.

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data used in logistic regression models can be found in